This R code includes the analyses and figures of the SMART manuscript.
It includes analyses of Nora Liebers, Peter-Martin Bruch, Junyan Lu, Miguel Hernandez and Tobias Terzer.
1 Setup
library(knitr)
options(digits=3, width=80)
opts_chunk$set(echo=TRUE,tidy=FALSE,include=TRUE)
1.1 Loading libraries
library(ggpmisc)
suppressMessages(library(RColorBrewer))
suppressMessages(library(ggrepel))
suppressMessages(library(gghighlight))
suppressMessages(library(ggbeeswarm))
suppressMessages(library(cowplot))
suppressMessages(library(gridExtra))
suppressMessages(library(drc))
suppressMessages(library(tidyverse))
detach("package:dplyr")
suppressMessages(library(dplyr))
library("survminer")
library("survival")
library("maxstat")
library("forestmodel")
library("forestplot")
library("patchwork")
library("table1")
library("glmnet")
library("corrplot")
library("writexl")
#maintain function of specific R packages
rename <- dplyr::rename
count <- dplyr::count
filter <- dplyr::filter
1.2 Define figure themes
fontsize=18
fontsize_small=10
starsize <- 5
statsize <- 5
theme_figures<- theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
axis.line = element_line(size=0.6),
axis.ticks = element_line(size=1.5),
legend.text = element_text(size=fontsize),
legend.title = element_text(size=fontsize+2))
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
theme_figures_wo_legend<- theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
axis.line = element_line(size=0.6),
axis.ticks = element_line(size=1.5),
legend.position = "none")
theme_figures_x_angle_45<-theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize, angle=45, vjust=1, hjust=1),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
axis.line = element_line(size=0.6),
axis.ticks = element_line(size=1.5),
legend.text = element_text(size=fontsize),
legend.title = element_text(size=fontsize+2))
theme_figures_x_angle_45_wo_legend<-theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize, angle=45, vjust=1, hjust=1),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
axis.line = element_line(size=0.6),
axis.ticks = element_line(size=1.5),
legend.position = "none")
theme_figures_x_angle<-theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize, angle=90, vjust=0.5, hjust=1),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
axis.line = element_line(size=0.6),
axis.ticks = element_line(size=1.5),
legend.text = element_text(size=fontsize),
legend.title = element_text(size=fontsize+2))
theme_figures_x_angle_wo_legend<-theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize, angle=90, vjust=0.5, hjust=1),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
axis.line = element_line(size=0.6),
axis.ticks = element_line(size=1.5),
legend.position = "none")
# used colors in combination with alpha=0.85
## primary colors
red1<-"orangered3"
gray1<-"slategray4"
blue1<-"deepskyblue4"
## secondary colors
red2<-"orangered4"
blue2<-"deepskyblue3"
1.3 Define functions
1.3.1 Functions for Table 1
my.render.cont <- function(x) {
with(stats.apply.rounding(stats.default(x), digits=3), c("",
"Median [Min, Max]"=sprintf("%s [%s, %s]", MEDIAN, MIN, MAX)))
}
my.render.cont_IQR <- function(x) {
with(stats.apply.rounding(stats.default(x), digits=3), c("",
"Median [Min, Max]"=sprintf("%s [%s, %s]", MEDIAN, IQR)))
}
my.render.cat <- function(x) {
c("", sapply(stats.default(x), function(y) with(
y,sprintf("%d (%0.0f %%)", FREQ, PCT))))
}
1.3.2 Functions for Analyses
#standard or Welch's two-sample t-test
t_test <- function(yourTab, var_eq=T){
pVal <- t.test(metric ~ response, yourTab, var.equal = var_eq)$p.value
}
m_reg <- function(yourTab){
#pvalue from the response comparison after accounting for tumor infiltration
return( summary( lm(AUC ~ response + t_infiltration, data = yourTab) )$coefficients[2,4] )
}
###Dose-response curves
fit_log_curve <- function(df, drugname){
#Fits a 5-parameter logistic model for each in-vivo response and returns:
# 1) coefficients (b,c,d,e,f) for each model.
# 2) r^2 between predicted and real values for each model.
# 3) a df with 1k predicted viability values
# for each in-vivo response. This can then be used for building a smooth dose-response curve.
fun1 <- function(df, drugname){
# Fit individual curves (5 parameter logistic worked best)
# parameters and meanings
# b, c, d, e, f
# slope, lower, upper, ED50, asymmetry
mR <- drm(normVal ~ concentration, data = df[df$response == "R",], fct = LL2.5(),
type = "continuous", upperl = c(NA, 1, NA, 1, NA), lowerl = c(NA, 0, NA, 0, NA))
mPD <- drm(normVal ~ concentration, data = df[df$response == "PD",], fct = LL2.5(),
type = "continuous", upperl = c(NA, 1, NA, 1, NA), lowerl = c(NA, 0, NA, 0, NA))
#get param coefficients
coefR <- mR$coef
coefPD <- mPD$coef
#evaluate fit
r2_mR <- round( cor(predict(mR),df[df$response == "R",]$normVal, use = "pairwise.complete.obs")^2, 2)
r2_mPD <- round( cor(predict(mPD),df[df$response == "PD",]$normVal, use = "pairwise.complete.obs")^2, 2)
#predict 1k concentrations to smoothen curves in line plot (dose-response curve)
df_conc <- data.frame(concentration = seq(min(df$concentration), max(df$concentration), length.out=1000))
predi_R <- predict(mR, newdata = df_conc )
predi_PD <- predict(mPD, newdata = df_conc )
return( list( coefR,
coefPD,
r2_mR,
r2_mPD,
tibble(name = rep(drugname, 2000),
response = c(rep("R", 1000), rep("PD", 1000)) ,
concentration = c( seq(min(df$concentration), max(df$concentration), length.out=1000),
seq(min(df$concentration), max(df$concentration), length.out=1000) ) ,
viab_pred = c(predi_R,predi_PD)) ) )
}
#function if error
errorfun <- function(){
#if fitting error, return a df with NAs
return( list( c(NA,NA,NA,NA,NA),
c(NA,NA,NA,NA,NA),
NA,
NA,
tibble(name = drugname, response = NA, concentration = NA, viab_pred = NA) ) )}
#call functions
out_lst <- tryCatch( fun1(df, drugname),
error = function(e){T} )
#handle error and return
if ( typeof(out_lst) == typeof(T) ){return( errorfun() )} else{return( out_lst )}
}
#Junyan's
#To calculate the area under the curve based on the trapezoidal rule
calcAUC <- function(value, conc) {
#value is a vector of normalized viability values
#conc is a vector of concentrations.
valueConc <- data.frame(viab=value, conc=conc)
valueConc <- valueConc[order(valueConc$conc),]
areaTotal <- 0
for (i in seq(length(value)-1)) {
areaEach <- (valueConc$viab[i] + valueConc$viab[i+1])*log(valueConc$conc[i+1]/valueConc$conc[i])*0.5
areaTotal <- areaTotal + areaEach
}
areaTotal/log(valueConc[nrow(valueConc),]$conc/valueConc[1,]$conc)
}
#Junyan's miscelaneous custom functions: library(jyluMisc)
#it creates a pdf doc where it arranges the plots contained in a plotList on a grid
makepdf <- function(x, name, ncol = 3, nrow = 2, figNum =NULL, width =20, height = 12) {
require(gridExtra)
#figs per page
if (is.null(figNum)) figNum = ncol*nrow
if (length(x) == 0) return(NULL)
#create empty pdf doc with dims
pdf(name, width = width, height = height)
#iterate over len of plots in plotList by figNum (e.g., 8 by 8 = 1, 9, 17...)
for (i in seq(1, length(x), by = figNum)) {
#print(names(x)[i])
j <- min(i + figNum - 1, length(x))
#arrange plots in grid
do.call(grid.arrange, c(x[i:j], ncol = ncol, nrow = nrow))
}
dev.off() %>% invisible #close pdf
}
# Functions for elastic net regression #
repeated_glmnet <- function(X, Y, alpha = alpha, type.measure = "auc", nlambda, nfolds, ncv = 1000, family = "binomial"){
folds.list <- lapply(1:ncv, function(x){return(balanced_folds(Y, nfolds))})
lasso.models <- lapply(X = folds.list, FUN = function(folds){
return(cv.glmnet(X, Y, alpha = alpha, type.measure = type.measure, intercept = TRUE, standardize = FALSE, nlambda = nlambda, foldid = folds, family = family))
})
tmp <- lapply(lasso.models, FUN = function(x){
lambda.min <- x$lambda.min
index.min <- x$lambda == lambda.min
cvm <- x$cvm[index.min]
coef <- x$glmnet.fit$beta[,index.min]
names(coef) <- rownames(x$glmnet.fit$beta)
return(list(lambda = lambda.min, cvm = cvm, coef = coef))
})
coefficients <- matrix(unlist(lapply(tmp, FUN = function(x){
return(x$coef)})), ncol = length(tmp), byrow = FALSE)
rownames(coefficients) <- names(tmp[[1]]$coef)
cvm <- sapply(tmp, FUN = function(x){
return(x$cvm)
})
lambda <- sapply(tmp, FUN = function(x){
return(x$lambda)
})
return(list(coef = coefficients, cvm = cvm, lambda = lambda))
}
give_repeated_lasso_results <- function(lasso.models, type.measure = "auc"){
coef_mat <-lasso.models$coef
median_coef <- apply(coef_mat, MARGIN = 1, FUN = function(x){
return(ifelse(any(x != 0), median(x[x != 0]), 0))})
proportion_selected <- apply(coef_mat, MARGIN = 1, FUN = function(x){
return(sum(x != 0)/length(x))})
cvms <- lasso.models$cvm
median_auc <- sapply(1:nrow(coef_mat), FUN = function(cur_index){
relevant_cols <- which(coef_mat[cur_index,] != 0)
return(ifelse(length(relevant_cols) > 0, median(cvms[relevant_cols]), 0.5))})
tmp.frame <- data.frame(Covariate_label = rownames(coef_mat), medianOR = median_coef, Proportion_selected = proportion_selected, median_AUROC = median_auc)
rownames(tmp.frame) <- tmp.frame$Covariate_label
tmp.frame$medianOR <- case_when(tmp.frame$Covariate_label %in% c("Vidaza", "Azaguanin", chemotherapies) ~ exp(tmp.frame$medianOR/10), TRUE ~ exp(tmp.frame$medianOR))
tmp.frame <- tmp.frame[order(tmp.frame$Proportion_selected, decreasing = TRUE),]
rownames(tmp.frame)
#Compute proportions matrix
coef_mat <- coef_mat[rownames(tmp.frame),]
mat.tmp <-matrix(0, nrow = nrow(coef_mat), ncol = nrow(coef_mat))
for(i in 1:nrow(mat.tmp))
{
relevant_cols <- which(coef_mat[i,] != 0)
for(j in 1:nrow(mat.tmp))
{
if(length(relevant_cols) != 0){
mat.tmp[i,j] <- sum(coef_mat[j, relevant_cols] != 0) /length(relevant_cols)
}
else
{
mat.tmp[i,j] <- 0
}
}
}
rownames(mat.tmp) <- colnames(mat.tmp) <- rownames(coef_mat)
diag(mat.tmp) <- proportion_selected[rownames(mat.tmp)]
return(list(overview_frame = tmp.frame, proportions.mat = mat.tmp))
}
# function to select folds in cv in a balanced fashion
balanced.folds.new <- function (y, nfolds = min(min(table(y)), 10))
{
totals <- table(y)
fmax <- max(totals)
nfolds <- min(nfolds, fmax)
nfolds = max(nfolds, 2)
folds <- as.list(seq(nfolds))
yids <- split(seq(y), y)
bigmat <- matrix(NA, ceiling(fmax/nfolds) * nfolds, length(totals))
for (i in seq(totals)) {
# cat(i)
if (length(yids[[i]]) > 1) {
bigmat[seq(totals[i]), i] <- sample(yids[[i]])
}
if (length(yids[[i]]) == 1) {
bigmat[seq(totals[i]), i] <- yids[[i]]
}
}
smallmat <- matrix(bigmat, nrow = nfolds)
smallmat <- pamr:::permute.rows(t(smallmat))
res <- vector("list", nfolds)
for (j in 1:nfolds) {
jj <- !is.na(smallmat[, j])
res[[j]] <- smallmat[jj, j]
}
return(res)
}
balanced_folds <- function (class.column.factor, cross.outer){
sampleOfFolds <- get("balanced.folds.new")(class.column.factor, nfolds = cross.outer)
permutated.cut <- rep(0, length(class.column.factor))
for (sample in 1:cross.outer){
# cat(sample, "\n")
permutated.cut[sampleOfFolds[[sample]]] <- sample
}
return(permutated.cut)
}
#Function to format floats
formatNum <- function(i, limit = 0.01, digits =1, format="e") {
r <- sapply(i, function(n) {
if (n < limit) {
formatC(n, digits = digits, format = format)
} else {
format(n, digits = digits)
}
})
return(r)
}
# scale matrics
mscale <- function(x, center = TRUE, scale = TRUE, censor = NULL, useMad = FALSE){
if (scale & center) {
if (useMad) {
x.scaled <- apply(x, 1, function(y) (y-median(y,na.rm = T))/meanAD(y))
} else {
x.scaled <- apply(x, 1, function(y) (y-mean(y,na.rm=T))/sd(y,na.rm = T))
}
} else if (center & !scale) {
if (useMad) {
x.scaled <- apply(x, 1, function(y) (y-median(y,na.rm=T)))
} else {
x.scaled <- apply(x, 1, function(y) (y-mean(y,na.rm=T)))
}
} else if (!center & scale) {
if (useMad) {
x.scaled <- apply(x, 1, function(y) y/meanAD(y))
} else {
x.scaled <- apply(x, 1, function(y) y/sd(y,na.rm = T))
}
} else {
x.scaled <- t(x)
}
if (!is.null(censor)) {
if (length(censor) == 1) {
x.scaled[x.scaled > censor] <- censor
x.scaled[x.scaled < -censor] <- -censor
} else {
x.scaled[x.scaled > censor[2]] <- censor[2] #higher limit
x.scaled[x.scaled < censor[1]] <- censor[1] #lower limit
}
}
return(t(as.matrix(x.scaled)))
}
2 Data read in and preprocessing
2.1 File path
# Modify if necessary
filepath <-"../../"
2.2 Load processed data
load(paste0(filepath,"data/SMARTrial_data.RData"))
screenData<-left_join(screenData, df_drugs, by="Drug")
data_all<-left_join(data_all, df_drugs, by="Drug")
Explanation of all datasets:
- screenData: processed ex-vivo raw data
- data_all: includes in-vivo and ex-vivo data of all 80 patients after data quality control
- df_chemo: includes only patients with chemotherapy and evaluable response (n=46) after data quality control, includes AUC and normalized viability per concentration and drug
- df_drugs: list and metadata of drugs used ex-vivo
3 Overview drugs
3.1 Figure: overview drugs
# prepare data
drugs_overview <-df_drugs %>% dplyr::count(pathway, FDA_approved) %>%
left_join(., dplyr::count(df_drugs, pathway), by="pathway")
drugs_overview$pathway <-factor(drugs_overview$pathway, levels=c("ABL (BCR)", "ALK", "Apoptosis", "BCR", "Bromodomain", "CDK", "Cell cycle" , "Chemotherapy", "DDR", "FLT3", "Hedgehog", "HGF", "Histone acetyltransferase", "Histone deacetylase", "Histone methyltransferase", "IDH", "JAK/STAT", "MAPK", "Notch", "Nuclear export", "PI3K/AKT/mTOR", "PKC", "PLK", "Promiscuous" , "Proteasome" , "Stress response", "TLR", "TNF/NFKB" ),
labels=c("BCR/ABL", "ALK", "Apoptosis", "B cell receptor", "Bromodomain", "Cyclin dependent kinases", "Cell cycle control" , "Chemotherapy", "DNA damage response", "FLT3", "Hedgehog", "HGF", "Histone acetyltransferase", "Histone deacetylase", "Histone methyltransferase", "IDH", "JAK/STAT", "MAPK", "Notch", "Nuclear export", "PI3K/AKT/mTOR", "PKC", "PLK", "Other" , "Proteasome", "Stress response", "TLR", "TNF/NFKB" ))
# define figure theme
theme_tmp<-theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize, angle=45, vjust=1, hjust=1),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
axis.line = element_line(size=0.6),
axis.ticks = element_line(size=1.5),
legend.text = element_text(size=fontsize),
legend.title = element_text(size=fontsize+2),
legend.position = c(0.7, 0.9))
# plot
ggplot(data=drugs_overview, aes(x=reorder(pathway, -n.y), y=n.x, fill=as.character(FDA_approved)))+
geom_bar(stat = "identity", width=0.6)+
scale_fill_manual(values=c("#cbd1d5", "#f2a4a3"), name="Drug type", label=c("Clinical development/tool compound", "FDA approved"))+
ylab("\n\nNumber of drugs\n")+
xlab("")+
theme_tmp
ggsave(width = 12, height =8, path=(paste0(filepath,"docs")), filename = "Figure_1B_Compound_library.svg")
write_xlsx(rownames_to_column(drugs_overview, var="rowname"), path=paste0(filepath,"docs/Source data/","Fig1B",".xlsx"))
# remove(drugs_overview, theme_tmp)
3.2 Supplementary table: Drug list
For appendix: includes the pathway, target, Supplier of each drug
'df_drugs_supp<-df_drugs%>%select(Drug, pathway, target, Supplier)%>%filter(Drug!="Bortezomib + Cytarabine(800nM)")%>%rename(Compound="Drug")%>%rename(Main.targets_mode.of.action="target")%>%rename(Simplified.drug.class="pathway")
write_xlsx(df_drugs_supp,"docs/table_S1_drugs.xlsx")
remove(df_drugs_supp)'
## [1] "df_drugs_supp<-df_drugs%>%select(Drug, pathway, target, Supplier)%>%filter(Drug!=\"Bortezomib + Cytarabine(800nM)\")%>%rename(Compound=\"Drug\")%>%rename(Main.targets_mode.of.action=\"target\")%>%rename(Simplified.drug.class=\"pathway\")\n\nwrite_xlsx(df_drugs_supp,\"docs/table_S1_drugs.xlsx\")\nremove(df_drugs_supp)"
4 Patient characteristics
4.1 Table 1
# preparation
df_invivo$diagnosis_summarized <-factor(df_invivo$diagnosis, levels=c("AML", "ALL","DLBCL", "BL", "FL", "MCL", "CLL", "B-PLL", "T-PLL", "T-NHL" ),
labels=c("AML", "ALL", "DLBCL", "Burkitt lymphoma", "Follicular lymphoma", "Mantle cell lymphoma", "CLL", "B-PLL", "T-PLL", "T-cell lymphoma, NOS" ))
label(df_invivo$age) <-"Age (years)"
label(df_invivo$diagnosis_summarized) <-"Diagnosis"
label(df_invivo$time_from_ED_dich) <-"Time from initial diagnosis (months)"
label(df_invivo$lines_pretreatment_cat) <-"Prior lines of treatment"
label(df_invivo$tumor_infiltration) <-"Tumor infiltration of sample"
label(df_invivo$Treatment_type) <-"Prescribed treatment at study entry"
df_invivo$material <-factor(df_invivo$material, levels=c("PB", "BM", "LN"),
labels=c("Periphal blood", "Bone marrow", "Lymph node"))
label(df_invivo$material) <-"Tumor sample origin"
df_invivo<-df_invivo %>%
mutate(Treatment_type=gsub("mall molecules","mall molecule inhibitors",Treatment_type))
# table 1 printing
table1(~ Sex + age + diagnosis_summarized + time_from_ED_dich + lines_pretreatment_cat+ material +tumor_infiltration +Treatment_type, data=df_invivo, overall="Total", render.continuous=my.render.cont, render.categorical=my.render.cat)
4.1.1 Plot of diagnosis
Fig_2B<-df_invivo %>%
dplyr::count(diagnosis) %>%
mutate(diagnosis=factor(diagnosis, levels= arrange(dplyr::count(df_invivo, diagnosis),desc(n))$diagnosis )) %>%
ggplot(aes(x="", y=n, fill=diagnosis))+
geom_bar(stat="identity", width=1) +
coord_polar("y", start=0)+
geom_text(aes(x=1.25, angle=c(70,0,0,0,0,30,55,0,0,40), label = ifelse(n>1,paste0(diagnosis, ": ", n),"")), position = position_stack(vjust=0.5)) +
theme_void()+
theme(legend.position = "bottom")+
scale_fill_brewer(palette="Set3")+
labs(fill="Diagnosis")
write_xlsx(rownames_to_column(Fig_2B$data, var="rowname"), path=paste0(filepath,"docs/Source data/","Fig2B",".xlsx"))
ggsave(width = 6, height = 4, path=(paste0(filepath,"docs")), filename = "Figure_2_Diagnosis_Piechart.svg")
4.1.1 Plot of SMARTrial subcohorts
Fig2D<-tibble("Subcohort"=c("Chemotherapy cohort", "CLL cohort treated with \nVenetoclax or Ibrutinib"), "n"=c(46,18)) %>%
ggplot(aes(x="", y=n, fill=Subcohort))+
geom_bar(stat="identity", width=1) +
coord_polar("y", start=0)+
geom_text(aes(x=1, size=5, label = ifelse(n>1,paste0( n, " / 64"),"")), position = position_stack(vjust=0.5)) +
theme_void()+
theme(legend.position = "bottom", legend.text = element_text(size=12))+
scale_fill_manual(values=c(blue1, red1))+
guides(size="none", x="none")+
labs(fill="")
Fig2D
ggsave(width = 6, height = 4, path=(paste0(filepath,"docs")), filename = "Figure_2_Subcohort_Piechart.svg")
write_xlsx(rownames_to_column(Fig2D$data, var="rowname"), path=paste0(filepath,"docs/Source data/","Fig2D",".xlsx"))
4.2 Patient list
For appendix: includes clinical data for each patient
df_patients_supp<-df_invivo%>%select(patientID, diagnosis, age, Sex, material, tumor_infiltration, lines_pretreatment, treatment_spec, Response_def_parameter_final, Resp_protocol, Death, time_EFS)
df_patients_supp<-df_patients_supp%>%mutate(time_EFS=ifelse(Resp_protocol=="na", "na", time_EFS))
write_xlsx(df_patients_supp,path=paste0(filepath,"docs/table_S1_SMARTrial_patients.xlsx"))
remove(df_patients_supp)
4.4 Time from Diagnosis to Treatment initiation
my_colors <- c("#7fc97f","#beaed4", "#fdc086", "#ffff99", "#386cb0" , "#f0027f")
# add names so bar for 'c' gets fill, too
names(my_colors) <- c("ALL","AML", "CLL","DLBCL", "MCL", "T-PLL")
df_invivo %>%
mutate(time_report=as.numeric(time_report), time_diarel_treatment=as.numeric(time_diarel_treatment)) %>%
mutate(colorscale=ifelse(time_report>time_diarel_treatment, diagnosis, "Report in time for treatment")) %>%
ggplot(aes(x=time_report, y=time_diarel_treatment))+
geom_beeswarm(aes(color=colorscale), cex=1.7, size=3)+
# scale_color_brewer(palette = "Paired")+
# scale_color_manual(values=c("#7fc97f","#beaed4", "#fdc086", "#ffff99", "#386cb0" ,"#f0027f", "red"), limits=c("ALL","AML", "CLL","DLBCL", "MCL", "T-PLL"))+
scale_color_manual(values = my_colors,limits = c("ALL","AML", "CLL","DLBCL", "MCL", "T-PLL"), na.value = "black")+
theme_figures+
scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), limits = c(0,160), breaks = c(0,2,7,14,30,90))+
scale_x_continuous(trans=scales::pseudo_log_trans(base = 10), limits = c(0,20), breaks = c(0,2,7,14))+
theme(legend.position = "bottom")+
coord_fixed()+
geom_vline(xintercept = 7, linetype=2)+
labs(color="", x="Time from study inclusion to DRP report (days)", y="Time from diagnosis of relapse / occurence of \ntreatment indication to treatment initiation (days)" )
df_invivo %>%
filter(time_report>time_diarel_treatment) %>%
group_by(diagnosis) %>%
count()
## # A tibble: 6 × 2
## # Groups: diagnosis [6]
## diagnosis n
## <chr> <int>
## 1 ALL 1
## 2 AML 13
## 3 CLL 4
## 4 DLBCL 1
## 5 MCL 2
## 6 T-PLL 2
df_invivo %>%
mutate(time_report=as.numeric(time_report), time_diarel_treatment=as.numeric(time_diarel_treatment),
delaybyDRP=time_report-time_diarel_treatment) %>%
filter(delaybyDRP>0) %>%
group_by(delaybyDRP) %>%
count()
## # A tibble: 7 × 2
## # Groups: delaybyDRP [7]
## delaybyDRP n
## <dbl> <int>
## 1 1 6
## 2 2 3
## 3 4 6
## 4 5 5
## 5 6 1
## 6 7 1
## 7 8 1
df_invivo %>%
group_by(diagnosis) %>%
summarise(median.time.to.treatment=median(time_diarel_treatment), range=paste0(min(time_diarel_treatment), " - ", max(time_diarel_treatment)), count = n()) %>%
arrange(desc(count))
## # A tibble: 10 × 4
## diagnosis median.time.to.treatment range count
## <chr> <dbl> <chr> <int>
## 1 AML 2.5 0 - 15 34
## 2 CLL 21 0 - 144 25
## 3 MCL 8 0 - 49 5
## 4 DLBCL 14.5 2 - 27 4
## 5 T-PLL 2 1 - 6 4
## 6 FL 32 22 - 82 3
## 7 ALL 5 0 - 10 2
## 8 B-PLL 6 6 - 6 1
## 9 BL 8 8 - 8 1
## 10 T-NHL 22 22 - 22 1
4.5 Overview of the drug respone pattern
df_forClustering <- data_all %>% filter(!Resp_protocol %in% c("",NA,"na")) %>%
dplyr::rename(response = Resp_protocol, pretreatment = Pretreatment) %>%
mutate(frozen = screenData[match(patientID, screenData$patientID),]$frozen) %>%
mutate(material = ifelse(material == "BM","Bone marrow",
ifelse(material == "LN","Lymph node", "Peripheral blood")),
sample = ifelse(frozen == TRUE,"Frozen","Fresh"),
TP53 = ifelse(TP53 %in% "1","Mut", ifelse(TP53 %in% "0", "WT", NA)))
viabMat <- df_forClustering %>% distinct(patientID, Drug, AUC) %>%
pivot_wider(names_from = patientID, values_from = AUC) %>%
column_to_rownames("Drug") %>% as.matrix()
viabMat <- viabMat[,complete.cases(t(viabMat))]
viabMat.scaled <- mscale(viabMat, censor =4)
annoCol <- distinct(df_forClustering, patientID, diagnosis, response, pretreatment, material, sample, TP53) %>%
data.frame() %>% column_to_rownames("patientID")
#define annotaiton colors
colorAnno <- list(diagnosis = structure(RColorBrewer::brewer.pal(10,"Set3"),
names = c("AML", "ALL","DLBCL", "BL", "FL", "MCL", "CLL", "B-PLL", "T-PLL", "T-NHL" )),
TP53 = c(Mut = red2, WT = "white"),
pretreatment = c(Yes = red2, No = "white"),
sample = c(Fresh = "#386cb0", Frozen = "#f0027f"),
response = c(R = blue1, PD = red1, SD = gray1),
material = c(`Bone marrow` = "#7fc97f", `Lymph node` = "#beaed4", `Peripheral blood` = "#fdc086")
)
p<-pheatmap::pheatmap(viabMat.scaled, annotation_col = annoCol, scale = "row", clustering_method = "ward.D2",
color=colorRampPalette(c(blue2, "white", red2))(100),
show_colnames = FALSE, annotation_colors = colorAnno,
main = "Hierarchical clustering of ex vivo responses in SMARTrial cohort", silent = TRUE)
plot(p$gtable)
# ggsave("../../docs/Figure_S4.png", plot = p$gtable, height = 15, width =15)
ggsave("../../docs/Figure_ED3.pdf", plot = p$gtable, height = 15, width =15)
write_xlsx(rownames_to_column(as.data.frame(viabMat.scaled), var="rowname"), path=paste0(filepath,"docs/Source data/","ED3",".xlsx"))
5 Primary endpoint
5.1 Calculation of primary endpoint
df_invivo%>%summarize(
mean.timelabreport=mean(time_report),
median.timelabreport=median(time_report),
min.timelabreport=min(time_report),
max.timelabreport=max(time_report))
## # A tibble: 1 × 4
## mean.timelabreport median.timelabreport min.timelabreport max.timelabreport
## <drtn> <drtn> <drtn> <drtn>
## 1 3.98 days 3 days 2 days 17 days
#calculate rate with labreport during 7 days
df_prim_endpoint<-df_invivo%>% select(patientID, time_report) %>% dplyr::mutate(within7d=ifelse(time_report<=7, 1, 0)) %>% dplyr::mutate(time_to_analysis=ifelse(time_report<=7, "within 7 days", "more than 7 days")) %>%
dplyr::mutate(duration_analysis_cat=case_when(time_report<=3 ~ "within 3 days", (time_report>4 &time_report<=7) ~ "within 4-7 days", (time_report>=8 &time_report<=10) ~ "within 8-10 days", time_report>10 ~ "more than 10 days" ))
bintest <- binom.test(sum(df_prim_endpoint$within7d), nrow(df_prim_endpoint))
bintest
##
## Exact binomial test
##
## data: sum(df_prim_endpoint$within7d) and nrow(df_prim_endpoint)
## number of successes = 73, number of trials = 80, p-value = 6e-15
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.828 0.964
## sample estimates:
## probability of success
## 0.912
cat("For ", sum(df_prim_endpoint$within7d)," of ", length(df_prim_endpoint$within7d), "patients the primary endpoint was reached. Therefore, for ", (bintest$estimate*100),"% (",round(bintest$conf.int[1],3),",", round(bintest$conf.int[2],3),") of participants ex-vivo drug testing was performed within a week.")
## For 73 of 80 patients the primary endpoint was reached. Therefore, for 91.2 % ( 0.828 , 0.964 ) of participants ex-vivo drug testing was performed within a week.
remove(bintest)
#calculate rate with labreport during 3 days
df_invivo%>% select(patientID,time_report) %>% dplyr::mutate(within3d=ifelse(time_report<4, 1, 0)) %>% dplyr::count(within3d) %>% mutate(percent = n / sum(n)*100)
## # A tibble: 2 × 3
## within3d n percent
## <dbl> <int> <dbl>
## 1 0 30 37.5
## 2 1 50 62.5
5.2 Visualisation primary endpoint (barplot)
df_prim_endpoint<-df_prim_endpoint%>% dplyr::count(time_to_analysis) %>% mutate(percent = n / sum(n)*100)
df_prim_endpoint$time_to_analysis<-factor(df_prim_endpoint$time_to_analysis, levels=c("within 7 days", "more than 7 days"), labels=c("within 7 days", "more than\n7 days"))
Fig2C<-ggplot(data=df_prim_endpoint, aes(x=as.factor(time_to_analysis), y=percent, fill=time_to_analysis))+
geom_bar(stat="identity", width=0.5, alpha=0.85)+
xlab("")+
annotate("text", x=1, y=98, label= "73/80", size=6)+
annotate("text", x=2, y=15, label= "7/80", size=6)+
scale_y_continuous("Number of patients in %", limits = c(0, 100), breaks=c(0, 20, 40, 60, 80, 100), expand = expansion(mult = c(0, 0.05)))+
scale_fill_manual(values=c(red1, blue1), name="Time to report")+
theme_figures_wo_legend
ggsave(width = 6, height = 4, path=(paste0(filepath,"docs")), filename = "Figure_2B_Time_to_report.svg")
write_xlsx(rownames_to_column(Fig2C$data, var="rowname"), path=paste0(filepath,"docs/Source data/","Fig2C",".xlsx"))
# remove(df_prim_endpoint)
6 Pathways and drug difference between in-vivo response groups in chemotherapy patients
6.1 Prepare base data frame
df_chemo$Pathway <-factor(df_chemo$Pathway,
levels=c("ABL (BCR)", "ALK", "Apoptosis", "BCR", "Bromodomain", "CDK", "Cell cycle" ,
"Chemotherapy", "DDR", "FLT3", "Hedgehog", "HGF", "Histone acetyltransferase",
"Histone deacetylase", "Histone methyltransferase", "IDH", "JAK/STAT", "MAPK",
"Notch", "Nuclear export", "PI3K/AKT/mTOR", "PKC", "PLK", "Promiscuous" , "Proteasome",
"Stress response", "TLR", "TNF/NFKB" ),
labels=c("BCR/ABL", "ALK", "Apoptosis", "B cell receptor", "Bromodomain", "Cyclin dependent kinases",
"Cell cycle control" , "Chemotherapy", "DNA damage response", "FLT3", "Hedgehog", "HGF",
"Histone acetyltransferase", "Histone deacetylase", "Histone methyltransferase", "IDH",
"JAK/STAT", "MAPK", "Notch", "Nuclear export", "PI3K/AKT/mTOR", "PKC", "PLK", "Other" ,
"Proteasome", "Stress response", "TLR", "TNF/NFKB" ))
df_chemo_bckup <- rename(df_chemo, response = resp_protocol)
df_chemo <- filter(df_chemo, resp_protocol != "SD") %>%
#rename some cols
rename(response = resp_protocol) %>%
#AUC per pathway per patient
group_by(patientID, Pathway) %>% mutate(AUC_pathway_patient = mean(AUC)) %>%
ungroup()
6.4 Panel B: Significance in R vs PD comparisons by drug (t-tests).
6.4.1 Prepare data
plotTab <- select(df_p_drugAUC, Pathway, name, pval) %>% unique() %>%
mutate(log10_pval = -log10(pval), significant = ifelse(pval < 0.05,T,F),
significant = factor(significant, levels=unique(significant) ),
name = as.character(name))
plotTab[plotTab$name == "Abemaciclib (LY2835219)",]$name <- "Abemaciclib"
# drug_order<-plotTab%>%
# arrange(pval)%>%
# select(Pathway)%>%
# mutate(Pathway=as.character(Pathway))%>%
# unlist()%>%
# unique()
6.4.2 Plot
Sign_drugs_RvsPD<-ggplot(plotTab, aes(x=factor(Pathway, levels = Pathway_order), y=log10_pval, fill=significant)) +
geom_point(shape=21, color ="black", stroke=0.5, size=2.5, alpha=0.6, show.legend = F) +
scale_fill_manual(values = c("TRUE"=red1, "FALSE"="white")) +
geom_hline(yintercept = -log10(0.05), linetype="22", col="gray37") +
gghighlight(significant == T, use_direct_label = F,
unhighlighted_params = list(fill = NULL, color="gray") ) +
geom_text_repel(data = filter(plotTab, significant == T), aes(label=name),
size=fontsize_small-4, max.overlaps = Inf, min.segment.length = 0.4,
box.padding = 0.65, segment.colour="grey")+
ylab(expression(-log[10]*'('*italic(P)~value*')')) +
xlab("") +
theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize, angle=90, vjust=0.5, hjust=1),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
panel.grid.major.x = element_line(size = 0.5,linetype = "13"),
axis.ticks.x = element_line(size = 0.5, color = "black", linetype = "13"),
axis.ticks.length.x = unit(0.3, "cm"),
legend.position = "none")
## Warning: Tried to calculate with group_by(), but the calculation failed.
## Falling back to ungrouped filter operation...
## Warning: Could not calculate the predicate for layer 2; ignored
Sign_drugs_RvsPD
# remove(plotTab)
Genotype - phenotype associations
FLT3-ITD ratio & FLT3 inhibitors
# prepare data
tmp<-data_all%>%
# here all AML patients with an available information about FLT3-ITD mutation were considered, even if patients had no evaluable response they were included because here we have no correlation with in-vivo response
filter(diagnosis=="AML")%>%
filter(!is.na(FLT3_ITD_ratio), FLT3_ITD_ratio!="na")%>%
# only FLT3 inhibnitors are considered
filter(Drug=="Crenolanib" | Drug=="Gilteritinib" | Drug=="Quizartinib"| Drug=="Sorafenib"| Drug=="Midostaurin hydrate")%>%
mutate(Drug=ifelse(Drug=="Midostaurin hydrate", "Midostaurin\nhydrate", Drug)) %>%
select(patientID, Target, diagnosis, FLT3_ITD_ratio, Drug_response_mean_low, Drug)
# correlation plot
plot_S2A<-tmp%>%ggplot(aes(x=Drug, y=Drug_response_mean_low, color=FLT3_ITD_ratio))+
geom_beeswarm(size=3, alpha=1, cex=5) +
scale_color_gradientn(colors=c(blue1, red1), limits=c(0,1), name="FLT3-ITD ratio", breaks=c(0,0.5,1.0))+
coord_cartesian(ylim = c(0,1.3)) +
scale_y_continuous(name="Viability (AUC)")+
theme_figures+
stat_cor(aes(x=FLT3_ITD_ratio, y=Drug_response_mean_low), method="kendall", cor.coef.name = c("tau"), label.sep="\n", label.x.npc=0.4, label.y.npc = 0.1, size=5)+
facet_wrap(~ Drug, nrow=1, scales="free_x")+
theme(axis.text.x = element_blank(),
axis.ticks.x=element_blank(),
axis.title.x = element_blank())
plot_S2A
## Warning: The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
FLT3-TKD & FLT3 inhibitors
# prepare dataset
tmp<-data_all%>%
# here all AML patients with an available information about FLT3-ITD mutation were considered, even if patients had no evaluable response they were included because here we have no correlation with in-vivo response, only patients with FLT3-TKD mutation but not FLT3-ITD mutation are considered
filter(diagnosis=="AML")%>%
filter(!is.na(FLT3), FLT3!="na")%>%
filter(FLT3_TKD==1 | FLT3==0)%>%
# only FLT3 inhibnitors are considered
filter(Drug=="Crenolanib" | Drug=="Gilteritinib" | Drug=="Quizartinib"| Drug=="Sorafenib"| Drug=="Midostaurin hydrate")%>%
mutate(Drug=ifelse(Drug=="Midostaurin hydrate", "Midostaurin\nhydrate", Drug)) %>%
select(patientID, Target, diagnosis, FLT3, Drug_response_mean_low, Drug)
# plot
plot_S2B<-ggplot(tmp, aes(x=Drug,y=Drug_response_mean_low,color=FLT3))+
geom_beeswarm(size=2.5, alpha=1, cex = 1.8)+
scale_color_manual(values=c(blue1, red1), name="", label=c("wildtype", "FLT3-TKD\nmutated"))+
scale_x_discrete(name="")+
scale_y_continuous(name="Viability (AUC)", limits=c(0, 1.4), breaks=c(0, 0.25, 0.5, 0.75, 1, 1.25))+
geom_hline(aes(yintercept=1), colour="black",
linetype="dashed") +
theme_figures_x_angle_45+
guides(color=guide_legend(title="FLT3", nrow=1))+
theme(legend.position = "bottom", axis.title.x = element_blank() )
plot_S2B
remove(tmp)
IDH mutations and venetoclax
# prepare dataset
tmp<- data_all%>%
# here all AML patients with an available information about IDH mutation were considered, even if patients had no evaluable response they were included because here we have no correlation with in-vivo response
filter(diagnosis=="AML")%>%
filter(Drug=="Venetoclax")%>%
select(patientID, IDH1, IDH2, Drug, Pathway, AUC)%>%
mutate(IDH=case_when(IDH1==1 ~ "IDH1\nmutated", IDH2==1 ~ "IDH2\nmutated", IDH1==0 & IDH2==0 ~ "wildtype", IDH1=="na" & IDH2=="na"~ "NA" ))%>%
filter(IDH!="NA")
# plot
plot_S2C<-ggplot(tmp, aes(x=IDH,y=AUC,color=IDH))+
# geom_boxplot()+
geom_boxplot(data=summarise(group_by(tmp, IDH), mean=mean(AUC), sd=sd(AUC), .groups = "keep"),
aes(x=IDH,color=IDH,lower=mean-sd,upper=mean+sd,middle=mean,ymin=mean-sd,ymax=mean+sd),
stat="identity", inherit.aes=FALSE)+
geom_beeswarm(size=3, cex=3)+
scale_color_manual(values=c(red1, gray1, blue1))+
scale_y_continuous(name="Viability (AUC)", limits=c(0, 1.2), breaks=c(0, 0.25, 0.5, 0.75, 1))+
geom_hline(aes(yintercept=1), colour="black",
linetype="dashed") +
ggpubr::stat_compare_means(method = "t.test",paired=F, size=starsize, comparisons = list( c(1,3)))+
theme_figures_wo_legend+
xlab("")
plot_S2C
remove(tmp)
Nutlin-3a and TP53
# prepare dataset
tmp<-data_all%>%
filter(!is.na(TP53), TP53!="na")%>%group_by(patientID, TP53)%>%
filter(Drug=="Nutlin-3a")%>%mutate(TP53=case_when(TP53==1 ~ "TP53\nmutated", TP53==0 ~ "wildtype"))
# plot
plot_S2D<-ggplot(tmp, aes(x= TP53,
y=AUC,
color=TP53))+
geom_boxplot(data=summarise(group_by(tmp, TP53), mean=mean(AUC), sd=sd(AUC), .groups = "keep"),
aes(x=TP53,color=TP53,lower=mean-sd,upper=mean+sd,middle=mean,ymin=mean-sd,ymax=mean+sd),
stat="identity", inherit.aes=FALSE)+
geom_beeswarm(size=3, cex=3)+
scale_y_continuous(name="Viability (AUC)", limits=c(0, 1.2), breaks=c(0, 0.25, 0.5, 0.75, 1))+
geom_hline(aes(yintercept=1), colour="black",
linetype="dashed") +
scale_x_discrete(name="")+
scale_color_manual(values=c(red1, blue1))+
ggpubr::stat_compare_means(method = "t.test",paired=F, size=starsize, comparisons = list( c(1,2)))+
theme_figures_wo_legend
plot_S2D
remove(tmp)
Plot Figure
design1<-"
AC
BD
EE
FF
"
tp<-theme(plot.tag = element_text(size=22, face="bold"))
wrap_elements(plot_S2A+theme(legend.position = "bottom"))+tp+
wrap_elements( plot_S2B+theme(axis.text.x = element_text(angle=0, hjust=0.5, vjust=1)))+tp+
wrap_elements( plot_S2C)+tp+
wrap_elements( plot_S2D)+tp+
wrap_elements( waterfall_plot_Rvs_PD+theme(axis.text.x = element_text(angle=45, hjust=1, vjust=1, size=14)))+tp+
wrap_elements( Sign_drugs_RvsPD+theme(axis.text.x = element_text(angle=45, hjust=1, vjust=1, size=14)))+tp+
plot_layout(design=design1, heights = c(1,1,1,1), widths = c(1,0.5))+
plot_annotation( tag_levels = "A" )
## Warning: The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
ggsave(path=(paste0("../../docs/")), filename = "Figure_ED5.pdf", device = cairo_pdf, width=14, height=19.8)
## Warning: The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
write_xlsx(rownames_to_column(plot_S2A$data, var="rowname"), path=paste0(filepath,"docs/Source data/","ED5A",".xlsx"))
write_xlsx(rownames_to_column(plot_S2B$data, var="rowname"), path=paste0(filepath,"docs/Source data/","ED5B",".xlsx"))
write_xlsx(rownames_to_column(plot_S2C$data, var="rowname"), path=paste0(filepath,"docs/Source data/","ED5C",".xlsx"))
write_xlsx(rownames_to_column(plot_S2D$data, var="rowname"), path=paste0(filepath,"docs/Source data/","ED5D",".xlsx"))
write_xlsx(rownames_to_column(waterfall_plot_Rvs_PD$data, var="rowname"), path=paste0(filepath,"docs/Source data/","ED5E",".xlsx"))
write_xlsx(rownames_to_column(Sign_drugs_RvsPD$data, var="rowname"), path=paste0(filepath,"docs/Source data/","ED5F",".xlsx"))
6.5 Panel C: Dose-response curves and AUC box-plots for exemplary drugs.
Selected drugs
example_drugs <- c("Idarubicin", "Vindesine")
Fit curves
#data storage
df_smooth_curve <- data.frame()
for ( drugname in example_drugs ){
#prepare data
df <- filter(df_chemo, name == drugname ) %>%
select(patientID, response, name, concentration, normVal) %>%
arrange( desc(response) ) %>% mutate(response = factor(response, levels = unique(response)))
#call fun
res <- suppressWarnings(fit_log_curve(df, drugname))
#store data (1k predicted viability values for plotting a smooth curve)
df_smooth_curve <- rbind(df_smooth_curve, res[[5]] )
}
prepare data
plotTab <- select(df_chemo, Pathway, name, patientID, response, concentration, normVal) %>%
#PD points on top
arrange( desc(response) ) %>% mutate(response = factor(response, levels=unique(response)))
Fig_3A_Source<-list(plotTab, df_p_drugAUC)
names(Fig_3A_Source)<-c("Data Dose Resp","Data Boxplot")
write_xlsx(Fig_3A_Source, path=paste0(filepath,"docs/Source data/","Fig3A",".xlsx"))
make figs
plotList_DResp <- lapply( example_drugs ,function(drugname){
#plot values: real with dots, predicted with line
plotTabSub <- filter(plotTab, name == drugname)
ciTab <- group_by(plotTabSub, response, concentration) %>%
summarise(meanVal = mean(normVal), sdVal = sd(normVal), n=length(patientID)) %>%
mutate(se = sdVal/sqrt(n), ci = 1.96*se)
p <- ggplot(NULL,aes(color=response, fill=response)) +
geom_line(data = df_smooth_curve[df_smooth_curve$name == drugname,],
aes(x= concentration, y= viab_pred), size=1) +
geom_point(data = plotTabSub ,
aes(x= concentration, y=normVal),
shape = 21, size=3.5, alpha = 0.3, color="black") +
geom_errorbar(data = ciTab, aes(x=concentration, y=meanVal, ymax = meanVal + ci, ymin=meanVal-ci),
width=0.1, size=1) +
scale_color_manual(values = c("R" = blue1, "PD"=red2)) +
scale_fill_manual(values = c("R" = blue2, "PD"=red1)) +
coord_cartesian(ylim = c(0,1.5)) +
scale_x_continuous(trans='log10') +
labs(title = paste0("",drugname) )+
xlab(expression("Concentration"~(mu*M))) + ylab("Viability") +
guides(fill="none", color="none") +
theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
axis.line = element_blank(),
panel.border = element_rect(size=1.5, fill = NA),
axis.ticks = element_line(size=1.5),
legend.text = element_text(size=fontsize),
legend.title = element_text(size=fontsize+2))
p
})
## `summarise()` has grouped output by 'response'. You can override using the
## `.groups` argument.
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `summarise()` has grouped output by 'response'. You can override using the
## `.groups` argument.
plotList_box <- lapply( example_drugs ,function(drugname){
#get the df
df_plot <- filter(df_p_drugAUC, name == drugname)
pval <- formatC(unique(df_plot$pval), digits = 2)
pText <- bquote(italic("P")~"="~.(pval))
#make plot
p <- ggplot(df_plot, aes(x=response, y=AUC)) +
geom_boxplot(aes(fill = response), width = 0.8, size=0.8, alpha = 0.85,
color = "black", outlier.shape = NA) +
scale_fill_manual(values = c("PD"=red1,"R"=blue1) ) +
geom_beeswarm(aes(color = response),cex=3, dodge.width=0.8, shape=21,
fill="slategray4", size=4, alpha= 0.7) +
ggplot2::annotate(geom = "text", label = pText, x=2.0, y=1.3, parse = FALSE, hjust=0)+
scale_color_manual(values = c("black", "black") ) +
geom_hline(yintercept = 1, linetype="22", col="gray30", size=0.7) +
coord_cartesian(ylim = c(0,1.5)) +
labs(title= "" ) +
xlab(expression(paste(italic("in vivo")," response"))) +
ylab("Viability (AUC)") +
guides(color="none", fill="none") +
theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
axis.line = element_blank(),
panel.border = element_rect(size=1.5, fill = NA),
axis.ticks = element_line(size=1.5),
legend.text = element_text(size=fontsize),
legend.title = element_text(size=fontsize+2))
return( p )
})
#merge both plot lists
plotList <- c(rbind(plotList_DResp, plotList_box))
get legend
plot figs
# plot_3C_tmp<-plot_grid(plotlist = plotList, nrow=4, ncol=2, rel_widths = c(1, 0.8))
#
#
# plot_3C<-plot_grid(plot_3C_tmp, legend_3C, ncol = 1, rel_heights = c(1, .1))
# plot_3C
design<-"
A#B#C#D
EEEEEEE
"
plot_3C<-
plotList[[1]]+
plotList[[2]]+
plotList[[3]]+
plotList[[4]]+
# plotList[[5]]+
# plotList[[6]]+
# plotList[[7]]+
# plotList[[8]]+
legend_3C+
plot_layout(design = design, widths = c(1,0.01,0.98,0.01,1,0.01,0.98), heights = c(1,0.2))
plot_3C
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
#remove(df, df_smooth_curve, plotList, plotList_box, plotList_DResp, plotTab, res, legend_3C)
7 Elastic net regression model
# Set a seed for reproducibility
set.seed(20191120)
#Modify drug names for formula function
tmp <- df_drugs %>% mutate(Drug = as.character(Drug), Drug = case_when(Drug == "8-Azaguanin" ~ "Azaguanin", Drug == "Azacytidine (Vidaza)" ~ "Vidaza",Drug == "Decitabine (Dacogen)" ~ "Decitabine", TRUE ~ Drug))
chemotherapies <- (tmp %>% filter(pathway == "Chemotherapy", !(Drug %in% c("Cisplatin", "Cyclophasphamide"))))$Drug
lasso_analysis_set <- data_all %>%
# filter only chemopatients
filter(chemo_pat==TRUE)%>%
filter(!is.na(Treatment_type), Treatment_type!="na", Resp_protocol %in% c("R", "PD"))%>%
# filter only patients with an evaluable in vivo response
filter(Response_not_evaluable!="1")%>%
# filter patients only with both plates available
filter(Count_plates==2) %>%
#restrict to chemotherapies, modify drug names
mutate(Drug = as.character(Drug)) %>%
mutate(Drug = case_when(Drug == "8-Azaguanin" ~ "Azaguanin", Drug == "Azacytidine (Vidaza)" ~ "Vidaza",Drug == "Decitabine (Dacogen)" ~ "Decitabine", TRUE ~ Drug)) %>%
filter(Drug %in% chemotherapies) %>%
select(patientID, Drug, AUC, Resp_protocol)
#Transform data to wide format for regression analysis
lasso_analysis_set <- lasso_analysis_set %>% tidyr::spread("Drug", AUC)
formula <- as.formula(paste("Resp_protocol ~ ", paste(chemotherapies, collapse = " + "), sep = ""))
x.bestr <- model.matrix(formula, lasso_analysis_set)
y.bestr <- lasso_analysis_set$Resp_protocol
lasso.bestr.models <- repeated_glmnet(X = x.bestr, Y = y.bestr, alpha = 0.3, nlambda = 100, nfolds = 3)
list.results.bestr <- give_repeated_lasso_results(lasso.bestr.models, type.measure ="auc")
#
# print.xtable(xtable(list.results.bestr[[1]]), include.rownames=FALSE, table.placement = "ht", type = "html")
# write.csv(list.results.bestr[[1]], file = "elastic_net_results.csv")
knitr::kable(list.results.bestr[[1]][,-1])
tabforFig3D<-list.results.bestr[[1]][1:6,-1]
colnames(tabforFig3D)<-c("Median OR", "Selected Proportion", "Median AUROC")
corrplot(list.results.bestr[[2]], type="upper", col=brewer.pal(n=8, name="RdYlBu"))
tabforFig3D<-tabforFig3D %>%
mutate(`Median OR`=round(`Median OR`,2),
`Selected Proportion`=round(`Selected Proportion`,2),
`Median AUROC`=round(`Median AUROC`,2))
Fig3D<- ggplot(rownames_to_column(tabforFig3D, var = "Drug"))+
annotate(geom = "table", size=fontsize_small-4, x = 0,y = 0, label = list(cbind(rownames_to_column(tabforFig3D, var = "Drug"))))+
theme_void()
Fig3D
tabforFig3D
## Median OR Selected Proportion Median AUROC
## Vincristine 0.96 1.00 0.85
## Idarubicin 0.97 0.99 0.85
## Vindesine 0.97 0.98 0.85
## Mitoxantrone 0.97 0.97 0.85
## Cladribine 0.97 0.88 0.84
## Fludarabine 0.98 0.53 0.83
write_xlsx(tabforFig3D, path=paste0(filepath,"docs/Source data/","Fig3B",".xlsx"))
8 EFS analyses
8.1 Drugs which were relevant in elastic net model
tmp<-data_all%>%
# filter only chemopatients
filter(chemo_pat==TRUE)%>%
filter(!is.na(Treatment_type), Treatment_type!="na")%>%
# filter only patients with an evaluable in vivo response
filter(Response_not_evaluable!="1")%>%
# filter only relevant drugs from elastic net
filter(Pathway=="Chemotherapy")%>%filter(Drug=="Vindesine" | Drug=="Idarubicin" | Drug=="Mitoxantrone" |Drug=="Vincristine" |Drug=="Cladribine"|Drug=="Fludarabine")%>%
# filter patients only with both plates available (to be consistent with elastic net model)
filter(Count_plates==2)
tmp$Drug<- as.factor(tmp$Drug)
tmp<-tmp%>%select(patientID, Drug, AUC, Event, time_EFS)
# AUC is scaled such that a unit change of the regressor corresponds to 10% change in cell viability.
tmp<-tmp%>%mutate(AUC_percent=AUC*100)
tmp<-tmp%>%mutate(AUC_10percent=AUC*100/10)
# loop to calculate HR + 95%CI + p for each drug
result<-data.frame(matrix(nrow =1, ncol=5))
colnames(result)<- c("drug", "HR", "lower", "upper","p_value_cox")
for(i in levels(tmp$Drug)){
tmp2<-tmp%>%filter(Drug==i)
cox<-coxph(Surv(time_EFS/365, Event) ~ (AUC_10percent),
data=tmp2)
cox_sum<-summary(cox)
result[i, 1] <-i
result[i, 2] <-cox_sum[["coefficients"]][, "exp(coef)"]
result[i, 3] <-cox_sum[["conf.int"]][, "lower .95"]
result[i, 4] <-cox_sum[["conf.int"]][, "upper .95"]
result[i, 5] <-cox_sum[["coefficients"]][, "Pr(>|z|)"]
}
result<-result[-which(rownames(result)==1),]
tablefig3D<-result[-which(rownames(result)==1),]
# plot results
knitr::kable(tablefig3D)
remove(tmp, tmp2, cox, cox_sum)
Plotting of Forest plot
result$drug<-factor(result$drug, levels=rev(c("Vincristine", "Idarubicin","Vindesine", "Mitoxantrone", "Cladribine", "Fludarabine")))
plot_3E<-ggplot(data=result ,aes(x=drug, y=HR, ymin=lower, ymax=upper,
)) + geom_pointrange() +
# scale_x_discrete(limits=result$drug,breaks=rev(levels(result$drug)), name="")+
scale_y_log10() +
xlab("")+
scale_color_manual(values="black", name="")+
coord_flip() +
geom_hline(aes(yintercept=1), colour="black", size=1.5,
linetype="dashed", alpha=0.3) +
# annotate("text", x = 1:nrow(result)+0.2, y = result$HR+0.003,
# label = ifelse( result$p_value_cox<0.001, paste0("p<","0.001"),
# paste0("p=", round(result$p_value_cox,3) ) ), size=fontsize_small-4)+
geom_text(nudge_x =.3, label=paste0("p=", round(result$p_value_cox,3)), size=fontsize_small-5)+
theme(panel.grid.minor = element_blank(),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
legend.position = "none",
panel.background = element_rect(fill = "white", color="black"),
panel.grid.major = element_blank(),
axis.ticks.x = element_line(size=1.5),
axis.ticks.y = element_blank())+
ylab("Hazard Ratio")
plot_3E
# remove(result)
8.2 EFS-KM curves for selected drugs
8.2.1 Vindesine ~ EFS in chemogroup
x<-"Vindesine"
# prepare data
tmp<-data_all%>%
# filter only chemopatients
filter(chemo_pat==TRUE)%>%
filter(!is.na(Treatment_type), Treatment_type!="na")%>%
# filter only patients with an evaluable in vivo response
filter(Response_not_evaluable!="1")%>%
# filter only relevant drugs from elastic net
filter(Pathway=="Chemotherapy")%>%filter(Drug=="Vindesine" | Drug=="Idarubicin" | Drug=="Mitoxantrone" |Drug=="Vincristine" |Drug=="Cladribine"|Drug=="Fludarabine")%>%
# filter patients only with both plates available (to be consistent with elastic net model)
filter(Count_plates==2)%>%
filter(Drug==x)
# perform maximally selected log rank statistics
m_EFS<-maxstat.test(Surv(time_EFS/365, Event) ~ AUC, data=tmp,
smethod="LogRank", pmethod="Lau94", minprop = 0.2,
maxprop = 0.8)
m_EFS
##
## Maximally selected LogRank statistics using Lau94
##
## data: Surv(time_EFS/365, Event) by AUC
## M = 2, p-value = 0.1
## sample estimates:
## estimated cutpoint
## 0.766
# divide patients into two groups based on cutpoint
tmp$drug_sensitive<-tmp$AUC<=m_EFS$estimate
# perform survival analyses
fit_EFS <- survfit(Surv(time_EFS/365, Event) ~ drug_sensitive, data=tmp)
surv_median(fit_EFS)
## strata median lower upper
## 1 drug_sensitive=FALSE 0.0932 0.0521 0.392
## 2 drug_sensitive=TRUE 0.4849 0.2301 NA
survdiff(Surv(time_EFS/365, Event) ~ drug_sensitive, data=tmp)
## Call:
## survdiff(formula = Surv(time_EFS/365, Event) ~ drug_sensitive,
## data = tmp)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## drug_sensitive=FALSE 15 14 7.29 6.19 8.44
## drug_sensitive=TRUE 28 19 25.71 1.75 8.44
##
## Chisq= 8.4 on 1 degrees of freedom, p= 0.004
# annotate p-value
pval <- filter(result, drug == "Vindesine")$p_value_cox
pval <- formatNum(pval, digits = 1)
pAnno <- bquote(italic("P")~"="~.(pval))
# plot KM curve
ggsurvplot_EFS_Vind<-ggsurvplot(fit_EFS, conf.int =FALSE, risk.table = TRUE, pval = FALSE,
title="Vindesine",
xlab ="Time (in years)",
xlim = c(0, 2.5),
break.time.by=0.5,
ylab ="\nEvent-free survival probability\n",
surv.median.line = "hv",
tables.theme=clean_theme(),
legend="none",
legend.title=(""),
legend.labs=c("less sensitive", "sensitive"),
palette = c(blue1, red1),
ggtheme=theme_figures)
plot_EFS_Vind <- ggsurvplot_EFS_Vind$plot + ggplot2::annotate("text",label=pAnno, x = 2, y=0.8, hjust =0, size=5)
plot_EFS_Vind
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
remove(tmp, x, m_EFS, fit_EFS, ggsurvplot_EFS_Vind)
8.2.2 Vincristine ~ EFS in chemogroup
x<-"Vincristine"
# prepare data
tmp<-data_all%>%
# filter only chemopatients
filter(chemo_pat==TRUE)%>%
filter(!is.na(Treatment_type), Treatment_type!="na")%>%
# filter only patients with an evaluable in vivo response
filter(Response_not_evaluable!="1")%>%
# filter only relevant drugs from elastic net
filter(Pathway=="Chemotherapy")%>%filter(Drug=="Vindesine" | Drug=="Idarubicin" | Drug=="Mitoxantrone" |Drug=="Vincristine" |Drug=="Cladribine"|Drug=="Fludarabine")%>%
# filter patients only with both plates available (to be consistent with elastic net model)
filter(Count_plates==2)%>%
filter(Drug==x)
# perform maximally selected log rank statistics
m_EFS<-maxstat.test(Surv(time_EFS/365, Event) ~ AUC, data=tmp,
smethod="LogRank", pmethod="Lau94", minprop = 0.2,
maxprop = 0.8)
m_EFS
##
## Maximally selected LogRank statistics using Lau94
##
## data: Surv(time_EFS/365, Event) by AUC
## M = 2, p-value = 0.3
## sample estimates:
## estimated cutpoint
## 0.722
# divide patients into two groups based on cutpoint
tmp$drug_sensitive<-tmp$AUC<=m_EFS$estimate
# perform survival analyses
fit_EFS <- survfit(Surv(time_EFS/365, Event) ~ drug_sensitive, data=tmp)
surv_median(fit_EFS)
## strata median lower upper
## 1 drug_sensitive=FALSE 0.0932 0.0767 0.603
## 2 drug_sensitive=TRUE 0.4849 0.2082 NA
survdiff(Surv(time_EFS/365, Event) ~ drug_sensitive, data=tmp)
## Call:
## survdiff(formula = Surv(time_EFS/365, Event) ~ drug_sensitive,
## data = tmp)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## drug_sensitive=FALSE 14 13 7.22 4.64 6.29
## drug_sensitive=TRUE 29 20 25.78 1.30 6.29
##
## Chisq= 6.3 on 1 degrees of freedom, p= 0.01
# annotate p-value
pval <- filter(result, drug == "Vincristine")$p_value_cox
pval <- formatNum(pval, digits = 1)
pAnno <- bquote(italic("P")~"="~.(pval))
# plot KM curve
ggsurvplot_EFS_Vinc<-ggsurvplot(fit_EFS, conf.int =FALSE, risk.table = TRUE, pval = FALSE,
title="Vincristine",
xlab ="Time (in years)",
xlim = c(0, 2.5),
break.time.by=0.5,
ylab ="\nEvent-free survival probability\n",
surv.median.line = "hv",
tables.theme=clean_theme(),
legend="none",
legend.title=(""),
legend.labs=c("less sensitive", "sensitive"),
palette = c( blue1, red1),
ggtheme=theme_figures)
plot_EFS_Vinc<-ggsurvplot_EFS_Vinc$plot + ggplot2::annotate("text",label=pAnno, x = 2, y=0.8, hjust =0, size=5)
plot_EFS_Vinc
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
remove(tmp, x, m_EFS, fit_EFS, ggsurvplot_EFS_Vinc)
8.2.3 Cladribine ~ EFS in chemogroup
x<-"Cladribine"
# prepare data
tmp<-data_all%>%
# filter only chemopatients
filter(chemo_pat==TRUE)%>%
filter(!is.na(Treatment_type), Treatment_type!="na")%>%
# filter only patients with an evaluable in vivo response
filter(Response_not_evaluable!="1")%>%
# filter only relevant drugs from elastic net
filter(Pathway=="Chemotherapy")%>%filter(Drug=="Vindesine" | Drug=="Idarubicin" | Drug=="Mitoxantrone" |Drug=="Vincristine" |Drug=="Cladribine"|Drug=="Fludarabine")%>%
# filter patients only with both plates available (to be consistent with elastic net model)
filter(Count_plates==2)%>%
filter(Drug==x)
# perform maximally selected log rank statistics
m_EFS<-maxstat.test(Surv(time_EFS/365, Event) ~ AUC, data=tmp,
smethod="LogRank", pmethod="Lau94", minprop = 0.2,
maxprop = 0.8)
m_EFS
##
## Maximally selected LogRank statistics using Lau94
##
## data: Surv(time_EFS/365, Event) by AUC
## M = 2, p-value = 0.5
## sample estimates:
## estimated cutpoint
## 0.548
# divide patients into two groups based on cutpoint
tmp$drug_sensitive<-tmp$AUC<=m_EFS$estimate
# perform survival analyses
fit_EFS <- survfit(Surv(time_EFS/365, Event) ~ drug_sensitive, data=tmp)
surv_median(fit_EFS)
## strata median lower upper
## 1 drug_sensitive=FALSE 0.132 0.0438 NA
## 2 drug_sensitive=TRUE 0.353 0.0986 1.66
survdiff(Surv(time_EFS/365, Event) ~ drug_sensitive, data=tmp)
## Call:
## survdiff(formula = Surv(time_EFS/365, Event) ~ drug_sensitive,
## data = tmp)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## drug_sensitive=FALSE 9 9 4.68 3.982 4.92
## drug_sensitive=TRUE 34 24 28.32 0.658 4.92
##
## Chisq= 4.9 on 1 degrees of freedom, p= 0.03
# plot KM curve
ggsurvplot_EFS_Clad<-ggsurvplot(fit_EFS, conf.int =FALSE, risk.table = TRUE, pval = FALSE,
title="Cladribine",
xlab ="Time (in years)",
xlim = c(0, 2.5),
break.time.by=0.5,
ylab ="\nEvent-free survival probability\n",
surv.median.line = "hv",
tables.theme=clean_theme(),
legend="none",
legend.title=(""),
legend.labs=c("less sensitive", "sensitive"),
palette = c( blue1, red1),
ggtheme=theme_figures)
ggsurvplot_EFS_Clad
plot_EFS_Clad<-ggsurvplot_EFS_Clad$plot
remove(tmp, x, m_EFS, fit_EFS, ggsurvplot_EFS_Clad)
9 Plot Figure 3
library("Cairo")
design<-"
AAAA
BBCC
DDDD
"
tp<-theme()
wrap_elements(plot_3C) +
wrap_elements(Fig3D) +
wrap_elements(plot_3E) +
wrap_elements(plot_3F) +
plot_layout(design=design, widths = c(1,1,.2,1.8), heights = c(1.2,.5,.7)) +
plot_annotation(title = "Figure 3", tag_levels = "A") & theme(title = element_text(size=30), plot.tag = element_text(size=22))
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
write_xlsx(as.data.frame(plot_3E$data), path=paste0(filepath,"docs/Source data/","Fig3C",".xlsx"))
ggsave(path=(paste0("../../docs/")), filename = "Figure_3.svg")
## Saving 16 x 16 in image
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
ggsave(path=(paste0("../../docs/")), filename = "Figure_3.pdf", device = cairo_pdf)
## Saving 16 x 16 in image
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
## Warning in is.na(x): is.na() auf Nicht-(Liste oder Vektor) des Typs 'language'
## angewendet
10 Burkitt lymphoma patient example (S005)
Burkitt Lymphoma: pretreatment with GMALL, during study treatment with R-DHAP –> PD with new liver lesions, then salvage treatment with MTX –> CR
# Prepare data: filter only for patient S005, label specific drugs which the patient had received to use these for colorcoding later in the plot
df_S005<-data_all%>%filter(patientID=="S005")%>%mutate(Drug_coding=case_when((Drug=="Doxorubicin" | Drug=="Cytarabine" ) ~ "refractory to treatment", Drug=="Pralatrexate" ~ "response to salvage treatment", (Drug!="Doxorubicin" &Drug!="Cytarabine" & Drug!="Pralatrexate") ~ "treatment not administered"))
# rename Drugs which includes the brand name in the original drug list
df_S005<-df_S005%>%mutate(Drug=ifelse(Drug=="Azacytidine (Vidaza)", "Azacytidine", Drug))
df_S005<-df_S005%>%mutate(Drug=ifelse(Drug=="Decitabine (Dacogen)", "Decitabine", Drug))
# define theme for figure
theme_figures_S005<- theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize, angle=45, vjust=1, hjust=1),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
axis.line = element_line(size=0.6),
axis.ticks = element_line(size=1.5),
legend.text = element_text(size=fontsize),
legend.title = element_text(size=fontsize+2),
legend.position ="bottom")
plot_S005_ESS<-df_S005%>%filter(Pathway=="Chemotherapy")%>%
ggplot(aes(x=reorder(Drug, ESS_AUC), y=ESS_AUC, fill=Drug_coding)) +
geom_bar(stat="identity", alpha=0.85)+
scale_x_discrete(name="")+
scale_y_continuous(name="Effect size score\n", limits=c(-0.5, 0.75))+
geom_hline(yintercept = 0)+
scale_fill_manual(name=expression(paste(italic("in vivo")," response")), values=c(blue1, red1, "gray84"))+
theme_figures_S005+
theme(axis.text.x =element_text(angle=90, vjust=0.25))+
guides(fill = guide_legend(nrow = 3))
plot_S005_ESS
ggsave( path=(paste0(filepath,"docs")), filename = "Figure_ED7A.svg")
## Saving 8 x 9 in image
write_xlsx(rownames_to_column(plot_S005_ESS$data, var="rowname"), path=paste0(filepath,"docs/Source data/","ED7A",".xlsx"))
remove(df_S005, plot_S005_ESS, theme_figures_S005)
11 CLL patients with venetoclax and ibrutinib treatment
11.1 Ibrutinib
theme_CLL_ibr<- theme_classic()+
theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
text=element_text(size=fontsize),
axis.text.y =element_text(size=fontsize),
axis.text.x =element_text(size=fontsize, angle=45, vjust=1, hjust=1),
axis.title.y =element_text(size=fontsize+2),
axis.title.x =element_text(size=fontsize+2),
axis.line = element_line(size=0.6),
axis.ticks = element_line(size=1.5),
legend.text = element_text(size=fontsize),
legend.title = element_text(size=fontsize+2),
legend.position = c(0.8, 0.2))
# prepare CLL ibrutinib dataset
df_CLL_ibr<-data_all%>%mutate(Ibrutinib=ifelse(grepl("Ibrutinib", treatment_spec), 1, 0))
df_CLL_ibr<-df_CLL_ibr%>%filter(Drug=="Ibrutinib")%>%filter(Ibrutinib==1)%>%filter(diagnosis=="CLL")%>%filter(Response_not_evaluable==0)
df_CLL_ibr$Resp_protocol <-factor(df_CLL_ibr$Resp_protocol, levels=c("R", "SD"),
labels=c("response","stable disease")) # patients had only R or SD, no patients with PD
plot_ibr<-df_CLL_ibr%>%ggplot(aes(x=reorder(patientID, -ES_AUC), y=-ES_AUC, fill=Resp_protocol)) +
geom_bar(stat="identity", alpha=0.75)+
scale_x_discrete(name="\nCLL patients")+
scale_y_continuous(name=expression(paste("Effect size (",italic("ex vivo"),")")), limits=c(-0.4, 0), breaks=c( -0.4, -0.2, 0), labels=c( 0.4, 0.2, 0))+
geom_hline(yintercept = 0)+
scale_fill_manual(values=c( blue1, red1), name=expression(paste(italic("in vivo")," response")))+
ggtitle("Ibrutinib")+
theme_CLL_ibr
# plot_ibr
remove(df_CLL_ibr, theme_CLL_ibr)
11.2 Venetoclax
# prepare CLL venetoclax dataset
df_CLL_ven<-data_all%>%mutate(Venetoclax=ifelse(grepl("Venetoclax", treatment_spec), 1, 0))
df_CLL_ven<-df_CLL_ven%>%filter(Drug=="Venetoclax")%>%filter(Venetoclax==1)%>%filter(diagnosis=="CLL")%>%filter(Response_not_evaluable==0)
df_CLL_ven$Resp_protocol<-factor(df_CLL_ven$Resp_protocol, levels=c("R"),
labels=c("response ")) # all patients had a R
plot_ven<-df_CLL_ven%>%ggplot(aes(x=reorder(patientID, -ES_AUC), y=-ES_AUC, fill=Resp_protocol)) +
geom_bar(stat="identity", alpha=0.75)+
scale_x_discrete(name="\nCLL patients")+
scale_y_continuous(name=expression(paste("Effect size (",italic("ex vivo"),")")), limits=c(-0.75, 0), breaks=c(-0.75, -0.5, -0.25, 0), labels=c(0.75, 0.5, 0.25, 0))+
geom_hline(yintercept = 0)+
scale_fill_manual(values=c(blue1), name=expression(paste(italic("in vivo")," response")))+
ggtitle("Venetoclax")+
theme_figures_x_angle_45_wo_legend
# plot_ven
remove(df_CLL_ven)
11.3 CLL VEN and Ibru cohort
design1<-"
AB
"
tp<-theme(plot.tag = element_text(size=22, face="bold"))
figure_CLL<-
wrap_elements(plot_ven)+tp+
wrap_elements(plot_ibr)+tp+
plot_layout(design=design1)
figure_CLL
ggsave(path=(paste0(filepath,"docs")), filename = "Figure_ED7C.svg")
## Saving 12 x 7 in image
# remove(tp, plot_ibr, plot_ven, figure_CLL)
write_xlsx(rownames_to_column(as.data.frame(plot_ven$data), var="rowname"), path=paste0(filepath,"docs/Source data/","ED7C",".xlsx"))
write_xlsx(rownames_to_column(as.data.frame(plot_ibr$data), var="rowname"), path=paste0(filepath,"docs/Source data/","ED7D",".xlsx"))
12 Quality assessment based on the type of sample material
(Analyses done by Miguel)
plate_mean_normVal of inner DMSO ctrls by sample_type: BM, LN, PB
anova_res <- filter(screenData, Drug== "DMSO", !edge) %>%
group_by(plateID, material) %>% summarise(plate_mean_normVal = mean(normVal)) %>% ungroup() %>%
mutate(material = factor(material, levels = unique(material)))
## `summarise()` has grouped output by 'plateID'. You can override using the
## `.groups` argument.
anova_res %>% DT::datatable(width = 700)
Show ANOVA res
m <- lm(plate_mean_normVal ~ material, data = anova_res)
summary(m)
##
## Call:
## lm(formula = plate_mean_normVal ~ material, data = anova_res)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16263 -0.01092 -0.00358 0.00641 0.27585
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.00255 0.00313 320.77 <2e-16 ***
## materialBM -0.00362 0.00763 -0.47 0.64
## materialLN 0.00192 0.00925 0.21 0.84
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0348 on 162 degrees of freedom
## Multiple R-squared: 0.00186, Adjusted R-squared: -0.0105
## F-statistic: 0.151 on 2 and 162 DF, p-value: 0.86
remove(anova_res, m)
Show plate median sdev and range
filter(screenData, Drug== "DMSO", !edge) %>% group_by(plateID, patientID) %>%
summarise(plate_sdev = round(sd(normVal),4) ) %>% ungroup() %>%
summarise(median_sdev = round(median(plate_sdev),4 ),
range_sdev = paste0( round(range(plate_sdev)[1],4) ,
" - ", range(plate_sdev)[2] ))
## `summarise()` has grouped output by 'plateID'. You can override using the
## `.groups` argument.
## # A tibble: 1 × 2
## median_sdev range_sdev
## <dbl> <chr>
## 1 0.0759 0.0293 - 0.6159
prepare data
plotTab <- select(screenData, plateID, plate_sdev, plate, material, frozen) %>% unique() %>%
#change names of categories
mutate(material = ifelse(material == "BM","Bone marrow",
ifelse(material == "LN","Lymph node", "Peripheral blood")),
frozen = ifelse(frozen == TRUE,"Frozen","Fresh"),
plate = ifelse(plate== "Acute","Myeloid","Lymphoid"))
show fig
p1<-ggplot(plotTab, aes(x = material, y = plate_sdev, fill = frozen, shape = plate)) +
geom_point(aes(color= plate), alpha = 0, size=2.5) +
geom_beeswarm(cex=1.5, size=2.5, alpha= 0.85, show.legend = F) +
scale_shape_manual(values = c("Myeloid" = 25, "Lymphoid" = 21)) +
scale_color_manual(values = c("Fresh"=red2,"Frozen"=blue1)) +
scale_fill_manual(values = c("Fresh"=red2,"Frozen"=blue1)) +
geom_hline(yintercept = 0.3, linetype="22", col="gray30", size=0.7) +
labs(color= "Material", shape = "Drug panel") +
ylab("Standard deviation (SD)\n") +
xlab("\nTumor sample origin") +
guides(fill = "none",
color = guide_legend(override.aes = list(alpha = 1), reverse = T, nrow=2),
shape = guide_legend(override.aes = list(alpha = 1), nrow=2)) +
theme_figures+
theme(legend.position = "right")
p1
ED_Fig_1A<-p1
# ggsave( path=(paste0(filepath,"docs")), filename = "Figure_S1_Quality_assessment.pdf")
write_xlsx(rownames_to_column(plotTab, var="rowname"), path=paste0(filepath,"docs/Source data/","ED1A",".xlsx"))
remove(plotTab)
12.1 Drug Drug Correlation
cor.mat <-
select(data_all, patientID, Drug, AUC) %>%
pivot_wider(names_from = Drug, values_from = AUC)%>%
column_to_rownames("patientID") %>%
rstatix::cor_mat()
cor.mat <- column_to_rownames(cor.mat, var="rowname")
breaks <- c(seq(-1, 1, length.out = 101)) %>% `names<-`(
colorRampPalette(c("#003DA5", "#2055B0", "#406EBC", "#6086C7", "#809ED2", "#9FB6DD", "#BFCFE9", "#DFE7F4", "white", "#F4E0E7", "#E9C2CF", "#DEA3B6", "#D3849E", "#C76586", "#BC476E", "#B12855", "#A6093D"))(101))
ED_Fig_1B<- pheatmap::pheatmap(cor.mat,
breaks = breaks,
color= names(breaks),
border_color=NA,
fontsize = 7,
filename=paste0(filepath,"docs/Figure_S2.pdf"),
width = 14,
height = 14)
write_xlsx(rownames_to_column(cor.mat, var="rowname"), path=paste0(filepath,"docs/Source data/","ED1B",".xlsx"))
# ggsave( path=(paste0(filepath,"docs")), filename = "Figure_S2.pdf")
# unique(data_all$patientID) %>% length()
12.2 Raw viability between fresh and frozen samples
screenData %>%
select(patientID, material, frozen) %>%
unique() %>%
group_by(material, frozen) %>%
count()
## # A tibble: 5 × 3
## # Groups: material, frozen [5]
## material frozen n
## <fct> <lgl> <int>
## 1 PB FALSE 35
## 2 PB TRUE 26
## 3 BM FALSE 8
## 4 BM TRUE 5
## 5 LN TRUE 9
Per Plate
CAVE: This treats plates as replicates
p1<-screenData %>%
filter(Drug %in% c("DMSO","PBS"), !edge) %>%
group_by(patientID, frozen, material) %>%
summarise(value=mean(value), sdev=mean(plate_sdev), .groups="keep") %>%
mutate(material = ifelse(material == "BM","Bone marrow",
ifelse(material == "LN","Lymph node", "Peripheral blood")),
frozen = ifelse(frozen == TRUE,"Frozen","Fresh")) %>%
ggplot(aes(x=frozen, y=value))+
geom_boxplot()+
geom_beeswarm()+
stat_compare_means(method="t.test", comparisons = list(c(1,2)))+
facet_wrap(~material)+
scale_y_log10()+
labs(y="Luminescence count\nin ATP assay", x="Sample pretreatment")+
theme_figures
p2<-select(screenData, plateID, plate_sdev, plate, material, frozen) %>% unique() %>%
mutate(material = ifelse(material == "BM","Bone marrow",
ifelse(material == "LN","Lymph node", "Peripheral blood")),
frozen = ifelse(frozen == TRUE,"Frozen","Fresh")) %>%
ggplot(aes(x=frozen, y=plate_sdev))+
geom_boxplot()+
geom_beeswarm()+
stat_compare_means(method="t.test", comparisons = list(c(1,2)))+
facet_wrap(~material)+
labs(y="Standard deviation (SD)", x="Sample pretreatment")+
theme_figures
p1+p2
Per Patient sample
p1<-screenData %>%
filter(Drug %in% c("DMSO","PBS"), !edge) %>%
group_by(patientID, frozen, material) %>%
summarise(value=mean(value), sdev=mean(plate_sdev), .groups="keep") %>%
mutate(material = ifelse(material == "BM","Bone marrow",
ifelse(material == "LN","Lymph node", "Peripheral blood")),
frozen = ifelse(frozen == TRUE,"Frozen","Fresh")) %>%
ggplot(aes(x=frozen, y=value))+
geom_boxplot()+
geom_beeswarm()+
stat_compare_means(method="t.test", comparisons = list(c(1,2)))+
facet_wrap(~material)+
scale_y_log10()+
labs(y="Luminescence count\nin ATP assay", x="Sample pretreatment")+
theme_figures
p2<-select(screenData, patientID, plateID, plate_sdev, plate, material, frozen) %>%
group_by(patientID, frozen, material) %>%
summarise(plate_sdev=mean(plate_sdev), .groups="keep") %>%
# unique() %>%
mutate(material = ifelse(material == "BM","Bone marrow",
ifelse(material == "LN","Lymph node", "Peripheral blood")),
frozen = ifelse(frozen == TRUE,"Frozen","Fresh")) %>%
ggplot(aes(x=frozen, y=plate_sdev))+
geom_boxplot()+
geom_beeswarm(cex=2)+
stat_compare_means(method="t.test", comparisons = list(c(1,2)))+
facet_wrap(~material)+
labs(y="Standard deviation (SD)", x="Sample pretreatment")+
theme_figures
p1+p2
14 Figure S3: Boxplots for all significant drugs: R vs SD vs PD
prepare data
#list of drugs that passed 0.05 (R vs PD). Arrange by pathway alphabetical (as fig 3B) and
#then by smallest pval.
example_drugs <- filter(df_p_drugAUC, pval <= 0.05) %>% arrange( Pathway , pval ) %>%
select(name) %>% unique()
#get the df
plotTab <- filter(df_chemo_bckup, name %in% example_drugs$name) %>%
select(-concentration, -concIndex, -normVal) %>% unique() %>%
mutate(plotTxt = paste0(Pathway," - ", name),
box_order = ifelse(response == "PD",0,ifelse(response == "SD", 1,2)) ) %>%
arrange( box_order ) %>% mutate(response = factor(response, levels=unique(response)))
remove(df_p_drugAUC)
boxplots
#make plot list
plotList <- lapply( example_drugs$name ,function(drugname){
p <- ggplot(plotTab[plotTab$name==drugname,], aes(x=response, y=AUC)) +
geom_boxplot(aes(fill = response), width = 0.8, size=0.3, alpha = 0.85,
color = "black", outlier.shape = NA) +
scale_fill_manual(values = c("PD"=red1,"SD"=gray1,
"R"=blue1) ) +
geom_beeswarm(aes(color = response),cex=2, dodge.width=0.8, shape=21,
fill="slategray4", size=2, alpha= 0.7) +
scale_color_manual(values = c("black", "black", "black") ) +
geom_hline(yintercept = 1, linetype="22", col="gray30", size=0.7) +
coord_cartesian(ylim = c(0,1.5)) +
labs(title= drugname ) +
xlab(expression(paste(italic("in vivo")," response"))) +
ylab("Viability (AUC)") +
guides(color="none", fill="none") +
theme_classic()+
theme( axis.line = element_blank(),
panel.border = element_rect(size=1, fill = NA),
plot.title = element_text(face='bold', hjust = 0.5, size = 14))
return( p )
})
See examples
cowplot::plot_grid(plotlist = plotList[1:8], ncol=2)
save in pdf
makepdf(plotList, paste0(filepath,"docs/Figure_S4_boxplots.pdf"), ncol =2,nrow = 4,
figNum = 8, width =8, height=12)
remove(plotList, plotTab, example_drugs, df_chemo_bckup)
15 Effect of Tumor Infiltration
Prepare data
df_tinf <- filter(screenData, plate_sdev < 0.3, Drug != "DMSO",
Pathway == "Chemotherapy", chemo_pat == T) %>%
select(patientID, diagnosis, treatment_type, chemo_pat, resp_protocol,
t_infiltration, plateID, Drug, concentration, normVal) %>%
#one viability per drug conc per patient
#(average, for drugs in both plates and replicate plates)
group_by(patientID, Drug, concentration) %>% mutate(normVal = mean(normVal)) %>%
select(-plateID) %>% unique() %>%
#AUC drug pat
group_by(patientID, Drug) %>% mutate(AUC = calcAUC(normVal, concentration)) %>%
select(-concentration, -normVal) %>% unique() %>%
#one viability per pat: mean AUC pat
group_by(patientID) %>% mutate(AUC_pat = mean(AUC)) %>% select(-Drug, -AUC) %>% unique() %>%
ungroup() %>%
mutate(t_infiltration = as.numeric(t_infiltration)) %>% #group_by(resp_protocol) %>%
#correlation test
mutate(corr_resp = round(cor(t_infiltration, AUC_pat, method = c("pearson")),2),
corr_pval = round(cor.test(x = t_infiltration, y = AUC_pat, method = c("pearson"))$p.value,2),
txt = paste0("R = ",corr_resp, "\np = ", corr_pval))
Plot
New Plot not subsetting by response group
p1<-ggplot(df_tinf, aes(x=t_infiltration, y=AUC_pat)) +
geom_point(size=2.5, show.legend = T, alpha=0.85) +
geom_smooth(method = "lm", se=F, show.legend = F)+
annotate("text", x = 60, y = 0.9,
label = c(unique(df_tinf[,]$txt)) ,
color=c("black"), size=5) +
coord_cartesian(xlim = c(50,100), clip = "off") +
xlab("\nTumor cell infiltration (in %)") + ylab("Viability (AUC)\n") +
theme_figures+
theme(aspect.ratio = 1)
p1
## `geom_smooth()` using formula = 'y ~ x'
plotTab<-filter(screenData, plate_sdev < 0.3, Drug != "DMSO",
Pathway == "Chemotherapy", chemo_pat == T) %>%
select(patientID, diagnosis, treatment_type, chemo_pat, resp_protocol,
t_infiltration) %>%
group_by(patientID, diagnosis, treatment_type, chemo_pat, resp_protocol,
t_infiltration) %>%
summarize(.groups = "keep") %>%
filter(resp_protocol%in%c("R","PD"))
plotTabmeanboxplot<-plotTab %>%
group_by(resp_protocol) %>%
summarise(mean=mean(t_infiltration), sd=sd(t_infiltration))
p2<- ggplot(plotTab, aes(x=resp_protocol,color=resp_protocol,y=t_infiltration))+
geom_boxplot(data=plotTabmeanboxplot, aes(x=resp_protocol,color=resp_protocol,lower=mean-sd,upper=mean+sd,middle=mean,ymin=mean-sd,ymax=mean+sd),stat="identity", inherit.aes
=FALSE)+
# geom_boxplot()+
geom_beeswarm(aes())+
stat_compare_means(comparisons = list(c(1,2)))+
coord_cartesian(ylim = c(0,100)) +
theme_figures+
guides(color="none")+
ylab("\nTumor cell infiltration (in %)") +
xlab(expression(paste(italic("in vivo")," response"))) +
scale_color_manual(values = c("R" = blue1, "PD"=red2))
p2
## Warning in wilcox.test.default(c(88, 70, 59, 80, 88, 90, 74), c(83, 93, : kann
## bei Bindungen keinen exakten p-Wert Berechnen
#eliminate conc-specific rows/cols
df_p_drugAUC <- select(df_chemo,
-concentration, -concIndex, -normVal, -AUC_pathway_patient) %>% unique() %>%
#add tumor infiltration variable
mutate(t_infiltration = as.numeric(screenData[match(patientID, screenData$patientID),]$t_infiltration) ) %>%
#"AUC" to "metric" for functions
mutate(metric = AUC) %>%
#perform tests
group_by(name) %>% nest() %>%
#t-test
mutate(pval_t = map(data, t_test),
#multivariate regression
pval_mreg = map(data, m_reg)) %>%
unnest( c(data, pval_t, pval_mreg) ) %>% ungroup() %>%
#eliminate patient data
select(Pathway, name, pval_t, pval_mreg) %>% unique() %>%
#add annotations
mutate(label_col = ifelse(pval_t > 0.05 & pval_mreg > 0.05, "01both_not_sig",
ifelse( pval_t < 0.05 & pval_mreg < 0.05, "02both_sig",
ifelse(pval_t < 0.05 & pval_mreg > 0.05, "03sig_t_only",
ifelse( pval_t > 0.05 & pval_mreg < 0.05, "04sig_mreg_only","") ) ) )) %>%
#plotting order of points
arrange(label_col)
df_p_drugAUC %>%
mutate_if(is.numeric, formatC, digits=2) %>% DT::datatable(width = 700)
p3<-ggplot(df_p_drugAUC, aes(x=-log10(pval_t) , y=-log10(pval_mreg), fill = label_col, color = label_col )) +
geom_point(size=3, shape = 21, alpha = 0.8, show.legend = T) +
scale_fill_manual(values = c("01both_not_sig"="#E18727FF","02both_sig"="deepskyblue4", "03sig_t_only"="darkred"),
labels = c("01both_not_sig"="p>0.05 in both models",
"02both_sig"="p<0.05 in both models",
"03sig_t_only"="p<0.05 without accounting\nfor tumor infiltration")) +
scale_color_manual(values = c("01both_not_sig"="#E18727FF","02both_sig"="deepskyblue4", "03sig_t_only"="darkred"),
labels = c("01both_not_sig"="p>0.05 in both models",
"02both_sig"="p<0.05 in both models",
"03sig_t_only"="p<0.05 without accounting\nfor tumor infiltration")) +
geom_abline(intercept = 0, slope = 1, linetype="22", col="gray30", size=0.7) +
xlab(expression(-log[10]*'('*italic(P)*'), tumor infiltration not considered')) +
ylab(expression(-log[10]*'('*italic(P)*'), accounting for tumor infiltration')) +
labs(fill = "Statistical significance",
color = "Statistical significance") +
theme_figures+
theme(legend.position = "right", aspect.ratio = 1)+
guides(fill=guide_legend(nrow=3))
p3
tp<-theme(plot.tag = element_text(size=22, face="bold"))
design2<-"
AB
CC"
wrap_elements(p1)+tp+
wrap_elements(p2)+tp+
wrap_elements(p3)+tp+
plot_layout(design=design2, heights = c(1,1), widths = c(2,1))+
plot_annotation(tag_levels = "A")
## `geom_smooth()` using formula = 'y ~ x'
## Warning in wilcox.test.default(c(88, 70, 59, 80, 88, 90, 74), c(83, 93, : kann
## bei Bindungen keinen exakten p-Wert Berechnen
ggsave( path=(paste0(filepath,"docs/")), filename = "Figure_ED6.pdf")
## Saving 12 x 12 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning in wilcox.test.default(c(88, 70, 59, 80, 88, 90, 74), c(83, 93, : kann
## bei Bindungen keinen exakten p-Wert Berechnen
write_xlsx(rownames_to_column(p1$data, var="rowname"), path=paste0(filepath,"docs/Source data/","ED6A",".xlsx"))
write_xlsx(rownames_to_column(p2$data, var="rowname"), path=paste0(filepath,"docs/Source data/","ED6B",".xlsx"))
write_xlsx(rownames_to_column(p3$data, var="rowname"), path=paste0(filepath,"docs/Source data/","ED6C",".xlsx"))
16 Midostaurin response of patients who received Midostaurin in vivo
# data_all %>%
# filter(diagnosis=="AML") %>%
# filter(treatment_spec%in%c("DA-Induction", "DA-Induction + Midostaurin" )) %>%
# group_by(patientID, FLT3, treatment_spec) %>%
# summarise() %>%
# # arrange(desc(FLT3))%>%
# arrange(treatment_spec)
data_all %>%
filter(diagnosis=="AML") %>%
filter(treatment_spec%in%c("DA-Induction", "DA-Induction + Midostaurin" )) %>%
filter(Drug%in%c("Daunorubicin","Cytarabine", "Midostaurin hydrate")) %>%
ggplot(aes(x=treatment_spec, y=AUC))+
geom_boxplot()+
geom_beeswarm(aes(color=FLT3))+
facet_wrap(~Drug)+
stat_compare_means(method="t.test")+
theme_figures_x_angle_45+
theme(legend.position = "bottom")+
xlab("In vivo therapy")+
ylab("Viability (AUC)")
Session info
sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8
## [3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
## [5] LC_TIME=German_Germany.utf8
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] Cairo_1.6-0 writexl_1.4.2 corrplot_0.92 glmnet_4.1-7
## [5] Matrix_1.5-1 table1_1.4.3 patchwork_1.1.2 forestplot_3.1.1
## [9] abind_1.4-5 checkmate_2.1.0 forestmodel_0.6.2 maxstat_0.7-25
## [13] survival_3.4-0 survminer_0.4.9 ggpubr_0.6.0 dplyr_1.1.1
## [17] lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0 purrr_1.0.1
## [21] readr_2.1.4 tidyr_1.3.0 tibble_3.2.1 tidyverse_2.0.0
## [25] drc_3.0-1 MASS_7.3-58.1 gridExtra_2.3 cowplot_1.1.1
## [29] ggbeeswarm_0.7.1 gghighlight_0.4.0 ggrepel_0.9.3 RColorBrewer_1.1-3
## [33] ggpmisc_0.5.2 ggpp_0.5.2 ggplot2_3.4.2 knitr_1.42
## [37] BiocStyle_2.26.0
##
## loaded via a namespace (and not attached):
## [1] TH.data_1.1-1 colorspace_2.1-0 ggsignif_0.6.4
## [4] ellipsis_0.3.2 markdown_1.5 gridtext_0.1.5
## [7] ggtext_0.1.2 rstudioapi_0.14 farver_2.1.1
## [10] MatrixModels_0.5-1 DT_0.27 fansi_1.0.4
## [13] mvtnorm_1.1-3 xml2_1.3.3 codetools_0.2-18
## [16] splines_4.2.2 cachem_1.0.7 polynom_1.4-1
## [19] Formula_1.2-5 jsonlite_1.8.4 broom_1.0.4
## [22] km.ci_0.5-6 cluster_2.1.4 pheatmap_1.0.12
## [25] BiocManager_1.30.20 compiler_4.2.2 backports_1.4.1
## [28] fastmap_1.1.1 cli_3.6.1 htmltools_0.5.5
## [31] quantreg_5.94 tools_4.2.2 gtable_0.3.3
## [34] glue_1.6.2 Rcpp_1.0.10 carData_3.0-5
## [37] jquerylib_0.1.4 vctrs_0.6.1 nlme_3.1-160
## [40] svglite_2.1.1 crosstalk_1.2.0 iterators_1.0.14
## [43] exactRankTests_0.8-35 xfun_0.38 timechange_0.2.0
## [46] lifecycle_1.0.3 gtools_3.9.4 rstatix_0.7.2
## [49] zoo_1.8-11 scales_1.2.1 ragg_1.2.5
## [52] hms_1.1.3 sandwich_3.0-2 SparseM_1.81
## [55] yaml_2.3.7 KMsurv_0.1-5 yulab.utils_0.0.6
## [58] sass_0.4.5 stringi_1.7.12 highr_0.10
## [61] foreach_1.5.2 plotrix_3.8-2 shape_1.4.6
## [64] commonmark_1.9.0 systemfonts_1.0.4 rlang_1.1.0
## [67] pkgconfig_2.0.3 evaluate_0.20 lattice_0.20-45
## [70] htmlwidgets_1.6.2 labeling_0.4.2 tidyselect_1.2.0
## [73] magrittr_2.0.3 bookdown_0.33 R6_2.5.1
## [76] magick_2.7.4 generics_0.1.3 multcomp_1.4-23
## [79] mgcv_1.8-41 pillar_1.9.0 withr_2.5.0
## [82] pamr_1.56.1 crayon_1.5.2 car_3.1-2
## [85] survMisc_0.5.6 utf8_1.2.3 tzdb_0.3.0
## [88] rmarkdown_2.21 data.table_1.14.8 digest_0.6.31
## [91] xtable_1.8-4 gridGraphics_0.5-1 textshaping_0.3.6
## [94] munsell_0.5.0 ggplotify_0.1.1 beeswarm_0.4.0
## [97] vipor_0.4.5 bslib_0.4.2
